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INTRODUCTION 

Finite element analysis offers great promise for reducing the 
empiricism now used extensively in polymer processing design, since 
it is well suited for modeling the complicated boundary conditions 
and material properties encountered in industrial practice. Thic 
paper will present the formulation we have developed for use in such 
processing analyses, and results which illustrate its use in some 
typical situations. Our element formulations are constructed on the 
premise that momentum convection can be neglected (polymer melt 
flows typically have very low Reynolds' numbers), but that convec- 
tive heat transfer may be significant (high Peclet numbers). 
Nonisothermal effects are considered important in polymer process- 
ing, due in part to the significant heating which may occur due to 
viscous dissipation, and also to the very strong influence of 
temperature on fluid viscosity. This paper will not discuss prob- 
lems in which the temperature dependence of viscosity is considered. 
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although our code does have the capability for ouch problems. Here 
we will treat the flow as Newtonian, with the flow field being 
coupled to the heat transfer equation only through the viscous heat 
generation. 

MODEL FORMULATION 


Governing Equations. 

The velocity and temperature fields are governed by the conserva- 
tion equations for momentum and energy as shown below ( a list of 
nomenclature appears at the end of the paper) : 

p(Du/Dt) » -VP + pV^u 

p c(DT/Dt) ~ kV^T + Q 

In our work the velocity field is constrained to be incompressible 
(Vu - 0) and the heat generation Q is given by the viscous dissi- 
pation (Q == x:Vn). Boundary conditions on velocity, stress, temper- 
ature and heat flux are also present on various portions of the 
boundary, and the ease with which these conditions are handled by 
the finite element procedure constitutes one of its most attractive 
features . 
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Finite Element Equations. 


As is described in detail in several texts (e.g. ref. 1)^ the 
finite element equations may be developed from the governing 
equations and their boundary conditions by recasting the equations 
in a Galerkin weighted residual formulation, applying Green's Theo- 
rem to lower the order of second derivatives and introduce natural 
boundary conditions, and finally by expressing the problem vari- 
ables in terms of an interpolation among their values at various 
nodes which are located at finite elements within the problem 
domain. 

In our formulation the nodal unknowns are taken as the velocity 
components and the temperature, so that: 

(u,T) = N^a^ 

where the i subscript refers to nodal values and the N functions are 
interpolation polynomials computed by standardized subroutines for 
isoparametric elements. The finite element equations can be written 
in the form: 


Ca + Ka = f 


where 
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The components in these global matrleeq ar*#=» j ^ 

y vwai inai;rices are assembled from corre 

.^ponding element eubmatrloes given in the following list.. 


c 


p 


N dfl 
Q 


(fluid inertia) 


OT . /■- 

J NpcN dfi (thermal inertia) 


k 




(m^B)^A (m^B) dQ 


(compressibility penalty) 


k 


|J _ 



B^DB dQ 


(viscosity) 


J (V%)kVN dfi 

Q 


(thermal conduction) 
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]f = I Npc iFvn dfi 


(thermal convection) 


= I NhN dr 


(boundary convection) 


Nt* dr 


(applied tractions) 


f = /*NQ dfi 

•I 

Q 


(heat generation) 




f = - / NhT^ d r 


(ambient temperature) 


We have chosen to compute various time- step values for transient 
problems using the "theta method" described in Reference 1. This 
algorithm may be written as: 


l(C/At) + + t(-C/At) - (l-e)Kla^ = (1-9 )£„ + 
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Here the n and n+1 subscripts refer to values at the current and next 
time, and 0 is a parameter between zero and one which determines the 
weight to be given the next time step. The algorithm is uncondi- 
tionally stable for 6 >.5, and 6=r667 corresponds to a Galerkin 
treatment with a linear interpolation over the time increment. 

All of this formulation is conventional, and a more detailed 
discussion of its underlying theory and computer implementation can 
be found elsewhere [1,2]. Some specific items might be mentioned 
here, however: (1) A penalty formulation is used to enforce incom- 
pressibility, rather than the usual velocity-pressure approach. 
This requires the use of double-precision computer arithmetic and a 
selectively reduced order of numerical integration for sufficient 
accuracy, but offers some reduction in programming effort and elimi- 
nates the need to compute pressure as an additional nodal variable. 
(2) The formulation assumes a full coupling between the viscous and 
thermal terms, with the resulting storage and manipulation of zeroes 
as seen in the above matrix equations. This coupling is unnecessary 
and inefficient for the partially coupled example problems 
described below, ’in which the flow field could be solved separately 
and then used in a single heat transfer solution. However, the 
coupled formulation is more general and we expect the majority of 
our work will require it. (3) We have coded a capability for either 
conventional Galerkin or "optimal” upwinding for the handling of the 
thermal convection term. The upwinding formulation is very conven- 
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lent, as it requires simply a one-point evaluation at a suitably 
shifted Gauss point within each element [ 2 ], but this approach is 
controversial and evidence can be cited as to its failure in some 
instances [ 3 ] , 

EXAMPLE PROBLEMS 

Entry Flow. 

Figure 1 shows the streamlines for a 4:1 entry flow which we have 
previously reported in greater detail [4]. Here a grid of 100 
four-node linear elements were used to model the upper symmetric 
half of the plane capillary , and a fully-developed Poiseuille 
velocity distribution was imposed at the reservoir entry as a bound- 
ary condition. The streamlines are identical with published 
experimental and numerical results, although the grid used here was 
not intended to be fine enough to capture the weak recirculation 
which develops in the stagnant corner of the reservoir. 

The temperature contours for convectionless flow are shown in 
Figure 2, which shows a hot region at the entrance of the capillary 
due to the combination of higher viscous dissipation at this region 
and a greater distance from cool boundaries to which heat may be 
conducted. These isotherms are normalized on the maximum centerline 
temperature T = UV^/3k expected for Poiseuille flow in the capil- 
lary. 
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The importance of thermal convection in this problem is given 
approximately by the Peclet number Pe “ ULpc/k, where wo may take U 
to be the maximum (centerline) velocity in the capillary and L ao the 
capillary half-height. Figure 3 plots the variation of temperature 
along the centerline for various values of Pe, and it can be seen 
that the effect of increased thermal convection is to sweep the 
cooler upstream flow particles well into the capillary, with a 
resulting lowering of the temperatures overall and a shift down- 
stream of the hot spot near the throat. The relatively coarse grid 
used in this problem produced unstable Galerkin results for Peclet 
numbers greater than approximately ten, and so the higher degrees of 
thermal convection were computed using the upwinding formulation. 
Further tests with refined grids should be completed to assess the 
accuracy of the upwinded solutions, although the plots in Figure 3 
appear reasonable. 

Transient Couette Flow. 

We have found that a capability for dynamic solutions is very 
useful not only in explicitly transient problems, but for a variety 
of processing flows. Rather than discretizing the entire bubble in 
a film blowing process, for instance, we follow a small annular 
strip dynamically as it travels from the die to the frost line. The 
dynamic algorithm is also useful for developing steady solutions to 
nonlinear problems, and we are presently working at modeling the 
flow of reactive fluids in this manner. 
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The performance of the time~otepping algorithm will be illus- 
trated here by means of some one-dimensional problems in which a 
single row of linear elements is used. Figure 4 shows such a strip 
which is used to model a Courttw flow in which the upper plate is set 
into motion impulsively, and in which the velocity gradients in the 
horizontal direction are zero. It is necessary in ouch problems to 
consider both the mesh spacing and the time step carefully, in order 
that both are able to capture sharp gradients which occur near the 
moving wall. It is noted that the mesh becomes much more refined 
near the upper plate in order to capture the boundary layer which 
develops there. We have also used a time step which begins very small 
and then increases logarithmically as the flow develops. The veloc- 
ity histories at three positions in the flow field are plotted in 
figure 4 and compared with the theoretical solution. It is clear 
that the computed solutions are not exact, but that the correct 
steady values are obtained at long times; improvement in accuracy 
could likely be obtained by considering further refinements in the 
mesh and the choice of time steps. 

Figure 5 shows the velocity and temperature profiles for a simi- 
lar problem, different only in that here the upper 15% of the fluid 
is a layer having a viscosity ten times that of the remaining 85%. 
This simulation is aimed at modeling the flow which may occur in 
coextrusion. At early times, the flow has not developed sufficient- 
ly to involve the low-viscosity fluid away from the upper plate. The 
thermal dissipation in the upper region is initially intense due to 
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the veiry high shear rates in the boundary layer, and a severe thermal 
spike develops. At longer times the flow is fully developed through 
the low-viscosity fluid? the high-viscosity region experiences only 
little shear flow and appears almost as a rigid boundary, The 
temperature profile in the low-viscosity fluid then becomes p«ara- 
bolic as expected in simple Couette flow. 

Graetz Flow. 

The Graetz or forced convection problem is a commonly used trial 
problem in which channel f,U;'W suddenly encounters a heated portion 
of the boundary. In the absence of flow the temperature distrib- 
ution is based only on conductive heat transfer considerations, but 
with significant flow the cooler upstream particles are swept down- 
stream; a thermal boundary layer develops at the heated boundary 
which grows gradually toward the channel center and the contours or; 
constant temperature are swept downstream relative to the 
conduction-only case. 

Figure 6 shows the isotherms for Peclet Numbers of 1.3 and 130, 
where upwinding was used to compute the high- convection case. The 
boundary conditions included an imposed pressure at the entrance of 
the channel, and the heated boundary was located a sufficient 
distance from the entry to permit the development of a Poiseuille 
flow. In this example the viscous generation of heat was assumed 
negligible, so heat is added only at the boundary. 
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CONCLUSIONS 


Thi 0 paper haa described some features of a model whieh is being 
developed in our laboratory to simulate a number of diverse polymer 
melt processing operations. We hope that this inexpensive and easi" 
ly implemented model can provide a means by which the designer's 
intuition might be expanded. Continuing work is aimed at increasing 
the number of process situations amenable to this type of modeling; 
among these are included reactive flows, free surfaces, and wall 
slip. 


NOMENCLATURE 

B Matrix of interpolation function derivatives 
c specific heat 

D Fluid viscosity matrix 

D/Du Substantive derivative 

h Coefficient of convective boundary heat transfer 

k Coefficient of thermal conduction 

L Characteristic length 

m Identity vector (1,1, 1,0, 0,0) 

n Unit normal vector 

N Interpolation function 

p Pressure 

T Temperature 

t Ambient temperature 


u Velocity vector 

V Rate of deformation tensor 

Boundary on which thermal convection occurs 
Boundary on which heat flux q* is specified 
r Boundary on which traction t* is specified 
Boundary or. which velocity u* is specified 
At time increment 

0 time step factor 

X Penalty coefficient (usually lO^y) 

y Newtonian viscosity 

p Density 

T Deviatoric stress tensor 

Q Problem domain (volume/ area, etc.) 

V Laplacian operator 
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Figure 3 - centerline temperatures in 4 j X entry 
flow, at various Peclet numbers. 








